function [fval, sol, EvaluationNum, GenerationNum] = differential_evolution(Problem,Algorithm)
%INPUT
%	Problem.
%       Dimension: scalar, the dimension of the decision variables
%       LowerBound: 1 -by- d vector, the lower bound of the decision variables
%       UpperBound: 1 -by- d vector, the upper bound of the decision variables
%       ObjFunc: the function handle of the objective funtions,
%           y = ObjFunc(x); x: n -by- d matrix; y: n -by- 1 vector
%   Algorithm.
%       PopulationSize: scalar, the population size
%       F: the mutation factor
%       CR: crossover index
%       StopCriteria: [X,Y]
%           [1,Y]: reach maximum budget Y
%           [2,Y]: reach the maximal generation number Y
%OUTPUT
%   fval: the optimal value
%   sol: the optimal solution
%   EvaluationNum: the total number of evalutions used
%   
%%
% PROBLEM SETTING
d = Problem.Dimension;
LB = Problem.LowerBound;
UB = Problem.UpperBound;
ObjFunc = Problem.ObjFunc;
% ALGORITHM ASSIGNMENT
PopulationSize = Algorithm.PopulationSize;
StopCriteria = Algorithm.StopCriteria;
F = Algorithm.F;
CR = Algorithm.CR;

N = (StopCriteria(1)==1)*StopCriteria(2) + (StopCriteria(1)~=1)*1000000; % the maximal budget size

%% INITIALIZATION
GenerationNum = 1;  % the current generation number
population = rand([PopulationSize,d]).*repmat((UB-LB),[PopulationSize,1])+repmat(LB,[PopulationSize,1]);
fitness = ObjFunc(population);
population = [fitness,population];  % [fitnesses, solution]
EvaluationNum = PopulationSize;  % the budget used
%% OPTIMIZATION
while EvaluationNum<N && (StopCriteria(1,1)~=2 || GenerationNum<StopCriteria(1,2))
    GenerationNum = GenerationNum+1;
    % mutation
    MutationVector = zeros(PopulationSize,d);
    for i = 1:PopulationSize
        r1 = ceil(rand().*PopulationSize);
        r2 = ceil(rand().*PopulationSize);
        while r2==r1
            r2 = ceil(rand().*PopulationSize);
        end
        r3 = ceil(rand().*PopulationSize);
        while r3==r1 || r3==r2
            r3 = ceil(rand().*PopulationSize);
        end
        MutationVector(i,:) = population(r1,2:end)+F.*(population(r2,2:end)-population(r3,2:end));
    end
    % crossover
    TrialVector = zeros(PopulationSize,d);
    for i = 1:PopulationSize
        u1 = (rand([1,d])<=CR);
        TrialVector(i,:) = population(i,2:end).*(1-u1) + MutationVector(i,:).*u1;
        u2 = ceil(rand().*d);
        TrialVector(i,u2) = MutationVector(i,u2);
    end
    fitness = ObjFunc(TrialVector);
    EvaluationNum = EvaluationNum+PopulationSize;
    betterid = find(fitness<population(:,1));
    population(betterid,:) = [fitness(betterid,:), TrialVector(betterid,:)];
end
[fval,opt_id] = min(population(:,1));
sol = population(opt_id,2:end);
end

